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Abstract. We consider a highly idealized model for El 
Nino/Southern Oscillation (ENSO) variability, as introduced 
in an earlier paper. The model is governed by a delay differ- 
ential equation for sea surface temperature T in the Tropical 
Pacific, and it combines two key mechanisms that participate 
in ENSO dynamics: delayed negative feedback and seasonal 
forcing. We perform a theoretical and numerical study of 
the model in the three-dimensional space of its physically 
relevant parameters: propagation period r of oceanic waves 
across the Tropical Pacific, atmosphere-ocean coupling 
and strength of seasonal forcing h. Phase locking of model 
solutions to the periodic forcing is prevalent: the local max- 
ima and minima of the solutions tend to occur at the same po- 
sition within the seasonal cycle. Such phase locking is a key 
feature of the observed El Nino (warm) and La Nina (cold) 
events. The phasing of the extrema within the seasonal cy- 
cle depends sensitively on model parameters when forcing is 
weak. We also study co-existence of multiple solutions for 
fixed model parameters and describe the basins of attraction 
of the stable solutions in a one-dimensional space of constant 
initial model histories. 



Keywords: Delay differential equations, El Nino, Ex- 
treme events. Fractal boundaries. Parametric instability. 



1 Introduction and motivation 

1 . 1 Key ingredients of ENSO theory 

The El-Nino/Southern- Oscillation (ENSO) phenomenon is 
the most prominent signal of seasonal-to-interannual climate 
variability. Its crucial role in climate dynamics and its socio- 
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(1199 ll) : iDiaz and Markgrai 



An international ten-year (1985-1994) Tropical-Ocean- 
Global- Atmos^here^TOGA^rogram greatly improved the 
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ing (iNeelin et all Il994[ll998h . and prediction (iLatif et al 



1998), theoretical model- 
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19941) of exceptionally strong El Nino events. It has con- 
firmed, in particular, that ENSO's significance extends far 
beyond the Tropical Pacific, where its causes lie. 

An important conclusion of this program was that — in 
spite of the great complexity of the phenomenon and the dif- 
ferences between the spatio-temporal characteristics of any 
particular ENSO cycle and other cycles — the state of the 
Tropical Pacific's ocean-atmosphere system could be char- 
acterized, mainly, by either one of two highly anticorrelated 
scalar indices. These two indices are a sea surface tempera- 
ture (SST) index and the Southern Oscillation Index (SOI): 
they capture the East-West seesaw in SSTs and that in sea 
level pressures, respectiv ely; see, for instance. Fig. 1 of 
Saunders and Ghill J2OO1I) . 

A typical version of the SST index is the so-called 
Nino-3.4 index, which summarizes the mean "anomalies" 
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Fig. 1. Temporal evolution of the NIN03.4 index that summa- 
rizes sea surface temperature (SST) anomalies in the region between 
170°W-120°W and 5°S-5°N. The time series is centered and nor- 
malized, but the horizontal lines do not represent the standard devi- 
ations: instead, they have have ordinates 1.5 and -1; see also Fig. [3] 



— i.e., the monthly-mean deviations from the climatological 
"normal" — of the spatially averaged SSTs over the region 



(170°W-120°W, 5° S-5° N) (IHurrell and Trenberth . 



Reynolds and Smithl.ll994 



Trenberth 



i\mi. 



1999; 



The evolution of this index since 1900 is shown in Fig.[TJ 
it clearly exhibits some degree of regularity, on the one hand, 
as well as numerous features characteristic of a determinis- 
tically chaotic system, on the other. The regularity mani- 
fests itself as the rough superposition of two don iinant oscil- 
lations — quasi-bienni al and quasi-quadrennial (IJiang et al. 



1995 



Ghil et al. 



20021) — accompanied by a near- symmetry 
of the local maxima and minima (i.e. , of the positive and neg- 
ative peaks). The lack of regularity h as been associated with 



the presence of a "Deyirs sta ircase" (IJin et al. 



1994, 



1996; 



Tziperman et al. 



1994 



19951) and does not preclude the su 



2008c ). 



perposition of stochastic effects as well (iGhil et al. 

While this study mainly focuses on local extrema (max- 
ima and minima) in our ENSO model, one must recall that 
the major El Nifios of 1982-83 and 1997-98 (see Fig.[T]) are, 
in fact, genuine extremes, i.e. rare events of unusually large 
magnitude. These climatic extremes and the related hydro- 
climatological impacts are part of the motivation for study- 
ing ENSO in general and for this study in particular. At the 
moment, the observational record contains too few of these 
truly extreme events to allow studying them by the methods 
of classical, i.e. statistical extreme value theory. We hope, 
therefore that the modeling approach developed in this study 
might prove useful in obtaining relevant statistical data for 
better understanding ENSO-related extreme events. 

To simulate, understand and predict such complex phe- 
nomena one needs a full hierarchy of models, from "toy" 
via intermediate to fully coupled general circulation mod- 



els (GCMs) (iNeelin et (2/.LIl998l:lGhil and RobertsonLl2000h . 
We focus here on a "toy" model, which captures a quali- 
tative, conceptual picture of ENSO dynamics that includes 
a surprisingly broad range of features. This approach al- 
lows one to gain a rather comprehensive understanding of the 
model's, and maybe the phenomenon's, underlying mecha- 
nisms and their interplay, at the cost of not capturing a full 
spatio-temporal picture of ENSO evolution. 

We consider the following conceptual ingredients that 
play a determining role in the dynamics of the ENSO phe- 
nomenon: (i) the Bjerknes hypothesis, which suggests a pos- 
itive feedback as a mechanism for the growth of an inter- 
nal instability that could produce large positive anomalie s 
of SSTs in the eastern Tropical Pacific (IBjerknessl Il969|) : 
(ii) delayed oceanic wave adjustments, realized in the form 
of eastward Kelvin and westward Ross by waves, that com 



pensate for Bjerknes 's positive fe edback (ISuarez and Schopi , 



1998 ) ; and ( i ii) seasonal forcing (lBattistiLll988l:lChang al. 



1994 
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Tziperman et al 



1994, 



2000). A more detailed discus- 



Ghil et al 



(l2008bh and 



sion of these ingredients is given by 
references therein. 

The past 30 years of research have shown that ENSO dy- 
namics is governed, by and large, by the interplay of the 
above nonlinear mechanisms and that their simplest version 



tems (Saunders and Ghil. 2001; Ghil etal 


, 2008a) and de- 


lav differential eauations (DDE) (ISuarez and School. 1998; 


Battisti and Hirst. 


1989; Tziperman et al . 1994). DDE mod- 



els provide a convenient paradigm for explaining interannual 
ENSO variability and shed new light on its dynamical prop- 
erties. So far, though, DDE model studies of ENSO have 
been limited to linear stability analysis of steady-state solu- 
tions, which are not typical in forced systems; case studies 
of particular trajectories; or one-dimensional (1-D) scenarios 
of transition to chaos, where one varies a single parameter 
while the others are kept fixed. A major obstacle for the 
complete bifurcation and sensitivity analysis of DDE mod- 
els lies in the complex nature of DDEs, whose analytical and 
numerical treatment is considerably harder than that of their 
ordinary differential equation (ODE) counterparts. 

1.2 Part 1 results and their physical interpretation 



Ghil et al 



(l2008bh took several steps toward a comprehen- 



sive analysis, numerical as well as theoretical, of a DDE 
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Fig. 2. Maximum and period maps for a warm initial history, (/)(t) = 1. (a) Maximum map, M — M{k, r) at 6 = 1; (b) Maxi mum map, 
M M (b, t) 2LtK = 10; (c) Period map, P = P{k, r) at 6 = 1; (d) Period map, P = P{b, r) Sit k = 10. Reproduced from 
with kind permission of Copernicus Publications on behalf of the European Geosciences Union (EGU). 
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model relevant for ENSO phenomenology. In doing so, they 
also illustrated the complexity of the structures that arise in 
its phase-and-parameter space for even such a simple model 
of climate dynamics. Specifically, the authors formulated a 
highly idealized DDE model for ENSO variability and fo- 
cused on the analysis of model solutions in a broad three- 
dimensional (3-D) domain of its physically relevant param- 
eters. They showed that this model can reproduce many 
scenarios relevant to ENSO phenomenology, including pro- 
totypes of El Nino and La Nina events; spontaneous inter- 
decadal oscillations; and intraseasonal activity reminiscent 
of Madden- Julian oscillations or westerly wind bursts. 

This model was also able to provide a good justification 



SSTs and trade winds ( 


Philander. 


1990: Diaz and Markeraf. 


1992: Jiane et al. 


1995 


: Ghil et al. 2002). with the 2-3-vear 



period arising naturally as the correct multiple of the sum of 
the basin transit t imes of Ke l vin and Rossby waves. An im- 
portant finding of iGhil et all (l2008b|) was the existence of re- 



gions of stable and unstable solution behavior in the model's 



parameter space; these regions have a complex and possibly 
fractal distribution of solution properties. 

Figure [2] illustrates the model's sensitive dependence on 
parameters in a region that corresponds roughly to actual 
ENSO dynamics. The figure shows the behavior of the global 
maximum M and period P of model solutions as a func- 
tion of three parameters: the propagation period r of oceanic 
waves across the Tropical Pacific, the atmosphere-ocean cou- 
pling strength and the amplitude b of the seasonal forcing; 
for aperiodic solutions we set P = 0. Although the model is 
sensitive to each of these three parameters, sharp variations 
in M and P are mainly associated with changing the delay r, 
which is plotted on the ordinate in all four panels of the fig- 
ure. In other words, the global maximum, in panels (a) and 
(b), as well as the period, in panels (c) and (d), may change 
more than twofold in response to a slight variation of r. 

This sensitivity is an important qualitative conclusion 
since in reality the propagation times of Rossby and Kelvin 
waves are affected by numerous phenomena that are not re- 
lated directly to ENSO dynamics. Moreover, the instabilities 
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disappear and the dynamics of the system becomes purely 
periodic, with period one year, as soon as the atmosphere- 
ocean coupling K, vanishes or the delay r decreases below a 
critical value; see Figs.[2t, b. Finally, the boundary between 
the domains of stable and unstable model behavior is clearly 
visible in Fig. [21 in the lower-right part of panels (b) and (d). 

The region below and to the right of this boundary con- 
tains simple period-one solutions that change smoothly with 
the values of model parameters. The region above and to 
the left is characterized by sensitive dependence on parame- 
ters. The range of parameters that corresponds to present-day 
ENSO dynamics lies on the border between the model's sta- 
ble and unstable regions. Hence, if the dynamical phenom- 
ena found in the model have any relation to reality. Tropi- 
cal Pacific SSTs and other fields that are highly correlated 
with them, inside and outside the Tropics, can be expected 
to behave in an intrinsically unstable manner; they could, in 
particular, change quite drastically with global warming. 

There are basically two approaches to ENSO dynamics 
( Neelin et al. i ll 994 Il998l) . both of which may be useful in 
extending the re sults of Pa r t 1 abo ve. The model consid- 
ered here and in 



Ghil et al. 



(I2008bh explains the complexi- 
ties of ENSO dynamics by the interplay of two oscillators: 
an internal, highly nonlinear one, due to a delayed feedback, 
and a forced, seasonal one. Our model thus falls within the 
strongly nonlinear, deterministic approach. 

An alternative approach attempts to explain several fea- 
tureas of ENSO dynamics by the action of fast, "weather" 
noise on a linear or very weakly nonlinear "slow" sys- 
tem, composed mainl y on th e up per ocean near the equa - 



tor. 



Boulanger et al\ (120041) and iLengaigne et all (120041) . 



among others, provide a comprehensive discussion of how 
weather noise could be responsible for the complex dynam- 
ics of ENSO, a nd, in particular , wheth er wind bursts trigger 



Savnisch et al. 



{ 20061) explore t his possibil- 



Ghil and Robertso n ( 2000) al- 



El Nino events, 
ity in a conceptual toy model, 
ready discussed the arguments about a "stochastic paradigm" 
for ENSO, with linear or only mildly nonlinear dynamics be- 
ing affected decisively by weather noise, v^. a "deterministi- 
cally chaotic parad igm," with decisively nonlinear dynamics. 



Ghil et al. 



(I2008d) have recently illustrated a way of combin- 
ing these two paradigms to obtain richer and more complete 
insight into climate dynamics in general. 

The present paper continues the study initiated in Part 1 
and focuses (i) on the multiplicity of model solutions for 
the same parameter values; and (ii) on the behavior of lo- 




Month 



Fig. 3. Histogram of temporal location of (a) warm and (b) cold 
events for the Nino-3.4 index. The event thresholds are shown by 
the dashed horizontal lines in Fig. [T] Notice the preferential occur- 
rence of both warm and cold events during the boreal winter. 



cal extrema in these solutions. In particular, we investi- 
gate the distribution in time of the model solutions' max- 
ima and minima; these extrema are directly connected to the 
strength and timing of the corresponding El Nino (warm) or 
La Nina (cold) events. The current analytic theory of DDEs 
does not allow one to easily answer many practically relevant 
questions about the behavior of even such apparently simple 
equations as our Eq. ([T]) below. The present study combines, 
therefore, general theoretical results about the existence and 
continuous dependence of solutions on parameters with ex- 
tensive numerical investigations. 

The rest of the paper is organized as follows. In Section O 
we summarize the model formulation from Part 1 , recall ba- 
sic theoretical results concerning this model's solutions, and 
briefly review details of the numerical integration method. 
Section [5] reports on the phase locking of solutions to the 
periodic forcing, namely on the tendency for the solutions' 
maxima and minima to each occur within a fixed, small in- 
terval of the seasonal cycle. Existence of multiple solutions 
and the attractor basins of the stable solutions are studied in 
Sect.m In Sect. Owe investigate the behavior of the model's 
local extrema, considered as a discrete dynamical system. A 
discussion of these results in Sect. [6l concludes the paper. 
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2 Model and numerical integration method 

2. 1 Model formulation and parameters 

Following Part 1 , we consider a nonlinear DDE with additive, 
periodic forcing, 

h\t) = -a tanh [f^h{t - r)] + 6 cos(27rcjt), (1) 

where h'{t) = dh{t)/dt, t > 0, and the parameters a, r, tz^ 6, 
and (J are all real and positive. Equation ([T]) is a simpli- 
fied one-delay ver sion o f the two-delay model considered by 
Tziperman et al. (119941) : it includes two mechanisms essen- 
tial for ENSO variability: a delayed, negative feedback via 
the function tanh(A>:z), and periodic external forcing. As 
shown in Part 1, these two mechanisms suffice to generate 
very rich behavior that includes several important features of 
more detailed models and of observational data sets. 

The function h(t) in ([T]) represents the thermocline depth 
deviations from the annual mean in the eastern Tropical 
Pacific; accordingly, it can also be interpreted roughly as 
the regional SST, since a deeper thermocline corresponds 
to less upwelling of cold waters, and hence higher SST, 
and vice versa. The thermocline depth is affected by the 
wind-forced, eastward Kelvin and westward Rossby oceanic 
waves. The waves' delayed effects are modeled by the func- 
tion tanh [nh{t — t)]\ the delay r is due to the finite velocity 
of these waves and it corresponds roughly to their combined 
basin-transit time. 

The particular form of the delayed nonlinearity plays 
an important role i n the behavior of a DDE model. 



Munnich et al. 



(119911) provided a physical justification for 
the monotone, sigmoid nonlinearity we adopt here. The 
parameter which is the linear slope of tanh(/^2;) at 
the origin, reflects the strength of the atmosphere-ocean 
coupling. The forcing term represents the seasonal cycle in 
the trade winds, with the strongest winds occurring in boreal 
fan. 

The DDE model ([U is fully determined by its five param- 
eters: feedback delay r, atmosphere-ocean coupling strength 

feedback amplitude a, forcing frequency a;, and forcing 
amplitude h. By an appropriate rescaling of time t and de- 
pendent variable /i, we let = 1 and a = 1. The remaining 
three parameters — and h — may vary, reflecting dif- 
ferent physical conditions of ENSO evolution. We consider 
here the same parameter ranges as in Part 1 of this study: 
0<r<2yr, 0</^<oo, 0<6<oo. 



To completely specify model ([B we need to prescribe 
some initial ''hist ory," /.g. th e behavior of h{t) on the in- 
terval [— r, 0), c/ iHald (1 1977b . In the numerical experiments 
of Sect. [3] below we assume, as in Part 1, that h{t) = 1, 
— r < t < 0, i.e. we start with a warm year. But in Sect. |4] 
we turn to a systematic exploration of the effect of the initial 
histories on the number and stability of solutions. 



2.2 Main theoretical result 

Consider the Banach space X 
uous functions h : [— r, 0) 

h e X SiS 



C{[-T, 0),R) of contin- 
^ and define the norm for 



h\\=sup{\h{t)lte[-r,0)}., 



where 



Nussbaum. 



denotes the absolute value in R (IHale . 



1977 



19981) . For convenience, we reformulate the DDE 
initial- value problem (I VP) in its rescaled form: 

h\t) = - tanh [/^ h{t - r)] + 6 cos(27r t), t > 0, (2) 
h{t) = (^(t) for t e [-r, 0), (l){t) e X. (3) 



Ghil et al. (l2QQ8bl) proved the fo l lowin g result, which 



Hale and Verduyn Lunell (Il993k and references 



follows from 
therein. 



Proposition 1 (Existence, uniqueness, continuous depen- 
dence) For any fixed positive triplet (r, b), the IVP (O- 
(O has a unique solution h(t) on [0, oo). This solution 
depends continuously on the initial data (j){t), delay r and 
the right-hand side of ([S]), considered as a continuous map 
/ : [0, T) X X ^ RJor any finite T. 

From Proposition [T] it follows, in particular, that the sys- 
tem (O-dS) has a unique solution for all time, which depends 
continuously on the model parameters (r, h) for any finite 
time. This result implies that any discontinuity in the solution 
profile, as a function of the model parameters, indicates the 
existence of an unstable solution that separates the attractor 
basins of two stable solutions. Our numerical experiments 
suggest, furthermore, that all stable solutions of (|2l)-([3]) are 
bounded and have an infinite number of zeros. 

2.3 Numerical integration 

The results in this Part 2 of our study are based on nu- 
merical integration of the DDE ([2l)-([3]). We emphasize 
that there are important differences between the numerical 
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Fig. 4. Seasonal phase locking of local extrema for model trajec- 
tories: (a-c) with period P = 1; (d-f) with period P = 7; and 
(g-i) aperiodic. The model solutions in panels (a, g) are shown 
in the stationary regime, after a sufficiently long transient, and the 
time axis is shifted so as to start from t = 0. The parameter values 
for these solutions are (a) r = 0.5, k, = 11, b = 2; (d) r = 0.56, 
K= 11, b = 1.4; and (g) r = 0.47, k = 10, b = 1.0. The scat- 
terplots of the points {h{ti), h{ti + 1)) in panels (6, e, h) use the 
values i = 0, 1, ... , 500, which correspond to to = 2500 and the 
parameter settings in panels (a,d,g), respectively The phase lock- 
ing is illustrated in panels (c, /, i), which give the /i-value of the 
local extrema — maxima shown as red filled circles and minima as 
blue squares — as a function of their position within the seasonal 
cycle, cp — t(mod 1). 
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Fig. 5. Seasonal phase locking of local extrema: cumulative results. 
Histogram of the phase Lp of the local maxima (red bars) and minima 
(blue bars) of model solutions with n = 2.0 (top panel) and n — 
11.0 (bottom panel). Each panel uses 10 000 individual solutions 
with parameters < r < 2 and < 6 < 10; see also Fig. [6] 



details of dde .solver, as well as a brief overview of other 
available DDE solvers, are given in Appendix C of Part 1. 



integration of DDEs and ODEs, and that these differences 
require developing special software; often the problem- 
specific modification of such software also becomes 
necessary. We used here the Fortran 90/95 DDE solver 



dde.solver of Shampine and Thompson (120061) . available 
at ,http://www.radford.edu/^hompson/ffddes/, Technical 



3 Seasonal phase locking of extrema 

A distinctive feature of the extreme ENSO phases — i.e., 
of the El Nino and La Nina events — is their occurrence 
during a boreal winter. This phenomenon is illustrated in 
Fig. [51 which shows the histograms of the monthly positions 
of unusually warm and unusually cold events, based on the 
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Nino-3.4 index of Fig. \T\ In our classification, El Ninos 
(see panel a) are those for which NINOS. 4 > 1.5, while La 
Ninas (see panel b) have NINOS. 4 < — 1. This asymmetry in 
the classification is due to the fact that extreme warm events 
are more intense but fewer in number than the extreme cold 



events ( 


Hoerlins et al. 


, 1997; 


Burgers and Stephenson. 


1999; 


Sardeshmukh et al . 2000: Kondrashov et al . 20051). Clearlv. 



the extreme events, both warm and cold, tend to occur during 
boreal winter. 



In discussing extrema, we distinguish between local and 
global ones. Recall that for a function h{t) specified on the 
interval [ti , ^2] , its global maximum (minimum) is defined as 
the point t such that h{t) is above (below) all the other values 
on that interval: h{t) > h{s), respectively h{t) < h{s), for 
all s G [^1,^2]. A local maximum (minimum) is a point t 
such that the corresponding value h{t) is above (below) all 
the values in a vicinity of t; for a sufficiently smooth func- 
tion, the latter definitions are equivalent to 

(i) h\t) = 0; and (ii) h'\t) < or h'\t) > 0, 

where h" = {h'Y is the second derivative of h{t). 

In this paper, we work with numerical solutions of the 
DDE problem ©-([S]) that are available only on a finite time 
interval [0, t/]; in addition, we eliminate the initial transient 
interval [0,to). We thus consider the global and local ex- 
trema of our solutions only on the interval [to^tf]. The global 
extrema thus defined might differ in certain cases from their 
counterparts on the interval [0, 00), for which our DDE is 
formally defined. The difference will only be noticeable for 
very-long-periodic, highly fluctuating solutions that are rel- 
atively rare in our model. Hence, the reduced definitions of 
the global and local extrema do not affect the main conclu- 
sions of our analysis. 

In this section, we study the phase (p of the local maxima 
and minima of the model solutions that obey ©-([S]). The 
main result, as we shall see, is that the model's extrema occur 
exclusively within a particular season. 

We start with several examples that illustrate the analysis 
in the rest of the section. Figure (111 shows a piece of model 
solution h{t) for r = 0.5, = 11, and b = 2. This solution 
has period P = 1, as illustrated in panel (6), which shows the 
scatterplot of the pairs {h{ti) ^ h{ti -\- 1)) for i = 0, 1, . . . and 
ti-^i = ti-\-l. Given the 1 -periodic character of the solution, 
all the points {h{ti)^h{ti + 1)) coincide. The choice of the 
starting point to will only affect the position of a single point 
in the panel (not shown). 
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Fig. 6. Seasonal phase locking of global extrema: parameter de- 
pendence. The plots show the phase ip of the global maxima of 
solutions of Eq. ^fovK, = 2.0 (top panel) and k = 11.0 (bottom 
panel); same number of solutions and parameter range as in Fig. [5] 
The rectangle in the bottom panel highlights the region blown up in 
Fig. El 

For each time epoch t we define its position (p within the 
seasonal cycle as = t(mod 1); the origin of the seasonal 
cycle in the forcing is taken in October, when the trade winds 
are strongest. Panel (c) shows the values of the local maxima 
(filled circles) and minima (squares) of h{t) as a function 
of their position (p within the seasonal cycle. The six other 
panels in Fig. H] show the results of a similar analysis for a 
solution with period P = 7 (panels d — f) and an aperiodic 
one (panels g — i). 

In all the examples of Fig.lH most of the local maxima are 
located within the first half of the annual cycle, i.e. in boreal 
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Fig. 7. Reversal in the phase locking of the maxima. The plot shows 
the seasonal cycle position Lp of the global maximum for 250 000 
solutions of Eq. for a€ = 11.0; it represents a blow-up of the 
region marked by a rectangle in the lower panel of Fig.[6l 

winter, while the local minima lie within the second half, i.e. 
in boreal summer. Moreover, the global maximum, as well 
as local maxima with large amplitudes, are always located 
within the (^-interval (0.15, 0.4), while the global minimum, 
as well as large-amplitude local minima, are always located 
within the interval (0.7, 0.95). We found this characteristic 
property of the model to hold for most of its solutions. 

To verify this model property, we analyzed the positions of 
the local extrema for a large number of individual solutions 
of Eq. (O within the parameter region (0<r<2,0<6< 
10) and at several values of n. The representative results are 
summarized in Figs. [5] and [6l where we used 10 000 individ- 
ual solutions for each value of n. Figured shows histograms 
of positions of the local extrema within the seasonal cycle, 
while Fig. [6l plots the position of the global maximum as a 
function of the model parameters r and h. 

The phase locking of the extrema to the seasonal cycle 
is present for most combinations of the physically relevant 
model parameters. Moreover, the local maxima tend to oc- 
cur, depending on the value of r, at either = 0.23 or 
(f = 0.27, while the local minima occur at (f = 0.73 or 
(p = 0.77. We notice that the cosine-shaped seasonal forc- 
ing vanishes at (p = 0.25 and (p = 0.75; hence the local 
maxima occur in the vicinity of zero forcing when the lat- 
ter decreases, and the local mimina occur in the vicinity of 
zero forcing when the latter increases. The offset in the posi- 




Time, t 



Fig. 8. Multiple stable solutions. Twenty trajectories that corre- 
spond to different initial histories (/)(t) = 00 collapse, after a tran- 
sient, onto four stable solutions. Two of these solutions are distinct, 
and the other two can be obtained from the latter by a time shift. 
Model parameters are r = 0.5, k = 10, and 6=1; see also Fig.|9l 

tion of the extrema from the point where the external forcing 
vanishes seems to be independent of the model parameters. 

As the atmosphere-ocean coupling parameter k, increases, 
yet another type of sensitive dependence on parameters sets 
in. Namely at low values of the external forcing, b < 1.5, 
"reversals" in the location of the local extrema do occur, with 
maxima suddenly jumping to boreal summer and minima to 
boreal winter. In Fig. [71 we zoom into one such reversal re- 
gion, marked by a rectangle in Fig. [6l The dark and light 
blue colors that occupy most of the region indicate that the 
global maximum of a model solution occurs in the first half 
of the annual cycle, while the red-to-yellow colors that ap- 
pear around r = 0.75 indicate that, within this "island," the 
global maximum jumps to the annual cycle's second half. 

4 Multiple solutions, stable and unstable 

The analysis in the previous section was carried out, as in 
Part 1, for the model ©-(Is]) with a fixed initial history, 
(/)(t) = 1. In this section, we study the model's solutions 
for distinct, yet still constant histories (/)(t) = (j)o. 

Naturally, different initial history values (j)o may result in 
different model solutions. This is illustrated in Fig. [8] for 
the parameter values r = 0.5, k, = 10, and b = 1. To 
produce this figure we used 20 distinct initial histories, with 
constant values that are uniformly distributed between 0o = 
—2 and 0o = 2; hence, at time t = there exist 20 distinct 
solutions. As time passes, those solutions are attracted by a 
smaller number of stable solutions so that, by t = 15, there 
are only four distinct solutions left, all of which have period 
P = 2. We notice furthermore that two of the remaining 
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"0 2 4 6 8 10 12 14 



Time, t 

Fig. 9. Solution profiles for multiple constant histories (/)(t) = 00- 
The top panel corresponds to point A = (r = 0.4, k = 11, 6 = 2) 
in parameter space, where there exists a unique stable solution. The 
bottom panel corresponds to B = (r = 0.5, k — 10, 6 = 1), the 
same point as in Fig. [S] here there exist two stable solutions and 
their attractor basins are bounded by horizontal discontinuity lines 
in the solution profile. The solutions are shown after a sufficiently 
long transient, and the origin of the time is shifted to start from zero; 
color bars indicate solution values, here as well as in Fig.fTol 

four solutions can be obtained by shifting the other two by 
one unit of time. 

In general, it is readily seen that — if the system ©-(13]) 
has solution x{t) — then x{t -\- k) with any integer k is also 
a solution. Hence, if x{t) is a solution with integer period 
P = k, then there are /c — 1 other solutions obtained from 
x{t) by an integer time shift. We will focus on solutions that 
cannot be obtained from each other by such a shift. Thus, we 
call two solutions x{t) and y{t) distinct if x{t) ^ y{t -\- k) 
for any positive integer k ^ P. 

Next we concentrate on the attractor basins of the model's 
stable solutions. Figure [9] shows the model's solution pro- 
files, after a suitable transient, for —lO<(j)o< 10, at two 
points in the model's parameter space: point A = {r = 
0.4, K = 1,6 = 2) in the top panel, and point B = {r = 
0.5, = 10, 6 = 1) in the bottom panel. Model behavior at 
point B was illustrated in Fig. [8l At point A the model has 
a unique stable solution that attracts all transient solutions 
as time evolves, so that the solution profile becomes constant 
along any vertical line, i.e. at any t = to in this type of figure. 

The model has two distinct stable solutions at point B: the 




Time, t 



Fig. 10. Multiple stable solutions. Solution profile for (a, c) dif- 
ferent initial histories (f){t) = 00, and (6, d) the corresponding dis- 
tinct solutions. For visual convenience, the trajectories are shifted 
to have their global maxima at t = 0. Panels (a, 6): model behavior 
at point C = (r = 0.5, k = 11,6 = 1.7842), where there ex- 
ist 2 distinct solutions; and panels (c, d): model behavior at point 
D ^ {t = 1.4579, = 11, 6 = 4), where there exist 61 distinct 
solutions. 

boundary between their attractor basins, as plotted on the real 
line of initial-history values corresponds to points of dis- 
continuity in the solution profiles. These points line up into 
straight horizontal lines in Fig. O one can see 8 horizontal 
lines of discontinuity in the solution profiles and there would 
thus appear to be 9 attractor basins. These basins correspond, 
however, as shown in Fig.Hl to only two stable solutions that 
are distinct from each other. 

Recall from Sect. 12. 21 that our solutions lie in the infinite- 
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Fig. 11. Local maxima (red) and minima (blue) of model solutions 
as a function of delay r; the other parameter values are fixed at 
K = 11 and 6 = 2. Notice the aperiodic regimes between periodic 
windows of gradually increasing period. 



dimensional Banach space X = r, 0), M), and that 

the solutions with constant initial histories do not span this 
space. By using such a particularly simple type of initial his- 
tories, we are merely exploring a 1-D manifold of solutions, 
parametrized by the scalar in the full space X. The in- 
tersection of the boundary between the attractor basins of the 
two stable solutions with this 1-D manifold gives the 8 lines 
of discontinuity seen in the bottom panel of Fig. [51 

Proposition 1 also implies that a discontinuity in the so- 
lution profile at suggests that there exists an unstable so- 
lution starting from = 0o- Hence, the boundary that 
separates the two attractor basins from each other is formed 
by unstable model solutions. This boundary is a manifold of 
codimension one in X, and Figure [9]re veals merely the inter- 
section of this manifold with the 1-D manifold of solutions 
that have constant initial histories. The presence of 8 such 
intersections suggests, in turn, that the boundary between the 
two attractor basins is a highly curved, but still smooth man- 
ifold. It is known for finite-dimensional problems that such 
boundaries can becom e quite complex and possibly fractal 



(|Grebogi et al. 



1987). 



Figure [Tol shows two slightly more complex situations 
along the same lines, namely one with still only two distinct 
solutions, having both period P = 2, but a more intricate pat- 
tern of solution profiles (panels a, 6), and one with 61 distinct 
solutions, having all P = 10 (panels c, d). For visual conve- 




0.5 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59 

Delay, t 



Fig. 12. Distribution of local maxima as a function of delay r within 
the interval 0.5 < r < 0.59; the other parameters are as in Fig. El 



nience we shift all the solutions so that their global maxima 
occur at t = 0. 



5 Dynamics of local extrema 

We focus here on the dynamics of the local extrema in the 
model solutions. For each solution h{t) we consider the se- 
quence of its local extrema {e^} := {h{ti), i = 1,2,...}, 
where h'{ti) = 0. The local maxima {Mi, i = 1, 2, . . . } are 
characterized by the additional condition that h"{ti) < 0, 
while at the local minima {rrii, i = 1,2,...} one has 
h'\U) > 0. 

Figure [TT] shows the position of the local extrema as a 
function of delay < r < 2 for fixed k, = 11 and b = 2. 
The figure illustrates convincingly the increase in complexity 
of model solutions as the delay r increases. For small delay 
values, < r < 0.5, each solution is a periodic sine-like 
wave with period P = 1, which contains a single maximum 
and a single mimimum within each cycle. 

Within the interval 0.6 < r < 0.8 the solutions become 
more complex: the solution period here is P = 3, and each 
cycle has three local maxima and three local minima. In gen- 
eral, the time elapsed between a local maximum and the next 
is an integer number; this effect is caused by the seasonal 
forcing, and the same is true for local minima. Often, the 
recurrence interval for extrema of the same kind is just unity 
and the number of local maxima (or minima) coincides with 
the period P of a given solution. 
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The period in Fig.fTTIincreases by jumps of 2, from P = 1 
to P = 3 and so on, as P = 2k-\-l. The transitions from one 
odd-periodic dynamics to the next are associated each time 
with a region of aperiodic behavior; e.g. the one from P = 1 
to P = 3 occurs in the interval 0.51 < r < 0.59. Thus, 
as r increases, the number of local extrema becomes larger 
and each increase in the number of extrema is preceded by a 
region of aperiodic, presumably chaotic behavior. 

Figure [121 zooms in on the distribution of local maxima 
within the first aperiodic region of Fig. [TTl namely 0.51 < 
r < 0.59. In this region, the r-intervals of aperiodic behav- 
ior alternate with shorter periodic windows: in the former 
the local maxima are distributed continuously within an in- 
terval, while in the latter several distinct local maxima oc- 
cur within a comparable range of values. This distribution 
of the maxima resembles the behavior of chaotic dynamical 
systen is in discrete time -— e.^., period dou bling for smooth 



maps (IFeigenbauml 1 1 97 8 : 



Kadanofi . 



19831) — and suggests 



that the model's aperiodic dynamics is in fact chaotic. An 
even richer behavior — with multiple, overlapping cascades 
— seems to emerge for 0.545 < r. 



6 Concluding remarks 

In the present paper we continued our study of a periodically 
force d delay differential equation (DDE) introduced by Ghil 
et al J2008bl) : the DDE © serves as a toy model for ENSO 
variability. We studied the model solutions numerically in a 
broad 3-D domain of physically relevant parameters: oceanic 
wave delay r, ocean- atmosphere coupling strength and 
seasonal forcing amplitude h. In this Part 2 of our investiga- 
tion, we focussed on multiple model solutions as a function 
of initial histories, and on the dynamics of local extrema. 

We found that the system is characterized by phase lock- 
ing of the solutions' local extrema to the seasonal cycle 
(Figs. |4] and [5]): solution maxima — i.e., warm events (El 
Ninos) — tend to occur in boreal winter, while local min- 
ima — i.e., cold events (La Ninas) — tend to occur in bo- 
real summer. The former model feature is realistic, since ob- 
served warm events do occur by-and-large in boreal winter; 
in fact, this property is one of the main features of the ob- 
served El Nino event s, having even given rise to th e name 



of the phenomenon (jPhilandei 



Diaz and Markgrai . 



199 



1990; 



Glantz et al 



1991 



The phase locking of cold events in the model to boreal 
summer is not realistic, since La Ninas also tend to occur in 
boreal winter, rather than in phase opposition to the warm 
ones; see again Fig. [3l It is not clear at this point which 
one of the lacking features of our DDE model gives rise to 
this unrealistic phase opposition; it might be the lack of a 
positive feedback mechanism, pres e nt wit h a separate, dis- 
tinct delay in the iTziperman et al.\ (119941) model. On the 
other hand, even GCMs with many more detailed features 
ma y have their warm events in entirely the wrong season; 



see 



Ghil and RobertsonI (l200Qh for a review. 



At the same time, for small-to-intermediate seasonal forc- 
ing b, the position of the global maxima and minima depends 
sensitively on other parameter values: it may exhibit signif- 
icant jumps in response to vanishingly small changes in the 
parameter values (Fig. [6]). In particular, an interesting phe- 
nomenon of "phase reversal" of the global extrema may oc- 
cur, cf. Fig. [71 

We explored a 1-D manifold of solutions for a set of given, 
prescribed points P = (r, a^, 6) in the model's parameter 
space. Such a manifold was generated, for each P, by so- 
lutions with constant initial histories (j){t) = 

We found multiple solutions coexisting for physically rel- 
evant values of the model parameters; see Figs.[8l-[TQl Some 
of these solutions are generated by shifting a single solution 
in time, using integer multiples of the period of the forcing, 
taken here to be unity. We have often found a set of k solu- 
tions so obtained from a single solution of period P = k. 

Typically, each stable solution has a bounded, but infinite- 
dimensional attractor basin in the solution space X described 
in Sect. 12.21 This attractor basin is separated from that of 
another stable solution by a manifold of codimension one, 
which is generated by unstable solutions (see Proposition 
1 and the following remarks). The intersections of such a 
manifold with the 1-D manifold of solutions explored herein 
appear as the straight horizontal lines in the solution-profile 
panels of Figs.l9land[TQl 

In Part 1 , we found that the solution period generally in- 
creases with the oceanic wave delay r. Figures [TT] and 
fT2l here provide much more detailed information: the pe- 
riod P of model solutions increases in discrete jumps, like 
{P = 2/c + 1, /c = 0, 1, 2, . . . }, separated by narrow, appar- 
ently chaotic "windows" in r. This increase in P is associ- 
ated with the increase of the number of distinct local extrema, 
all of which tend to occur at the same position within the 
seasonal cycle. The distribution of the maxima in Fig.fT2lre- 
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sembles in fac t the behavior of chaotic dynamical s ystems in 



discrete time (IFeigenbaum . 



1978 



Kadanoff. 



19831) and sug- 



gests that the model's aperiodic dynamics is in fact chaotic. 

It is quite interesting that, for plausible values of the de- 
lay r, the periods lie roughly between 2 and 7 years, a range 



strons warm events (Philander. 


1990 


; Glantz et al. 




1991; 


Diaz and Markgraf 


19921 Neelin et al 


. 1998L The sensitive 



dependence of the period on the model's external parameters 
(r, K^h) is consistent with the irregularity of occurrence of 
strong El Ni nos, and can help explain the difficulty in pre- 
dicting them (iLatif al\ 1 19941: lohil and Jiangi 119981) . 
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